! This code is part of
! RRTM for GCM Applications - Parallel (RRTMGP)
!
! Eli Mlawer and Robert Pincus
! Andre Wehe and Jennifer Delamere
! email:  rrtmgp@aer.com
!
! Copyright 2015-,  Atmospheric and Environmental Research,
! Regents of the University of Colorado, Trustees of Columbia University.  All right reserved.
!
! Use and duplication is permitted under the terms of the
!    BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause
! -------------------------------------------------------------------------------------------------
!>
!> ## Numeric calculations for gas optics. Absorption and Rayleigh optical depths, Planck source functions.
!>
!>   - Interpolation coefficients are computed, then used in subsequent routines. 
!>   - All applications will call compute_tau_absorption(); 
!>     compute_tau_rayleigh() and/or compute_Planck_source() will be called depending on the 
!>     configuration of the k-distribution. 
!>   - The details of the interpolation scheme are not particaulrly important as long as arrays including 
!>     tables are passed consisently between kernels. 
!>
! -------------------------------------------------------------------------------------------------

module mo_gas_optics_rrtmgp_kernels
  use mo_rte_kind,      only : wp, wl
  use mo_rte_util_array,only : zero_array
  implicit none
  private
  public :: interpolation, compute_tau_absorption, compute_tau_rayleigh, compute_Planck_source
contains
  ! --------------------------------------------------------------------------------------
  !> Compute interpolation coefficients
  !> for calculations of major optical depths, minor optical depths, Rayleigh,
  !> and Planck fractions
  subroutine interpolation( &
                ncol,nlay,ngas,nflav,neta, npres, ntemp, &
                flavor,                                  &
                press_ref_log, temp_ref,press_ref_log_delta,    &
                temp_ref_min,temp_ref_delta,press_ref_trop_log, &
                vmr_ref,                                        &
                play,tlay,col_gas,                              &
                jtemp,fmajor,fminor,col_mix,tropo,jeta,jpress) bind(C, name="rrtmgp_interpolation")
    ! input dimensions
    integer,                            intent(in) :: ncol,nlay
      !! physical domain size
    integer,                            intent(in) :: ngas,nflav,neta,npres,ntemp
      !! k-distribution table dimensions 
    integer,     dimension(2,nflav),    intent(in) :: flavor
      !! index into vmr_ref of major gases for each flavor
    real(wp),    dimension(npres),      intent(in) :: press_ref_log
      !! log of pressure dimension in RRTMGP tables 
    real(wp),    dimension(ntemp),      intent(in) :: temp_ref
      !! temperature dimension in RRTMGP tables 
    real(wp),                           intent(in) :: press_ref_log_delta, &
                                                      temp_ref_min, temp_ref_delta, &
                                                      press_ref_trop_log
      !! constants related to RRTMGP tables
    real(wp),    dimension(2,0:ngas,ntemp), intent(in) :: vmr_ref
      !! reference volume mixing ratios used in compute "binary species parameter" eta

    ! inputs from profile or parent function
    real(wp),    dimension(ncol,nlay),        intent(in) :: play, tlay
      !! input pressure (Pa?) and temperature (K)
    real(wp),    dimension(ncol,nlay,0:ngas), intent(in) :: col_gas
      !! input column gas amount - molecules/cm^2 
    ! outputs
    integer,     dimension(ncol,nlay), intent(out) :: jtemp, jpress
      !! temperature and pressure interpolation indexes 
    logical(wl), dimension(ncol,nlay), intent(out) :: tropo
      !! use lower (or upper) atmosphere tables 
    integer,     dimension(2,    ncol,nlay,nflav), intent(out) :: jeta
      !! Index for binary species interpolation 
#if !defined(__INTEL_LLVM_COMPILER) && __INTEL_COMPILER >= 1910
    ! A performance-hitting workaround for the vectorization problem reported in
    ! https://github.com/earth-system-radiation/rte-rrtmgp/issues/159
    ! The known affected compilers are Intel Fortran Compiler Classic
    ! 2021.4, 2021.5 and 2022.1. We do not limit the workaround to these
    ! versions because it is not clear when the compiler bug will be fixed, see
    ! https://community.intel.com/t5/Intel-Fortran-Compiler/Compiler-vectorization-bug/m-p/1362591.
    ! We, however, limit the workaround to the Classic versions only since the
    ! problem is not confirmed for the Intel Fortran Compiler oneAPI (a.k.a
    ! 'ifx'), which does not mean there is none though.
    real(wp),    dimension(:,       :,   :,    :), intent(out) :: col_mix
#else
    real(wp),    dimension(2,    ncol,nlay,nflav), intent(out) :: col_mix
      !! combination of major species's column amounts (first index is strat/trop)
#endif
    real(wp),    dimension(2,2,2,ncol,nlay,nflav), intent(out) :: fmajor
      !! Interpolation weights in pressure, eta, strat/trop 
    real(wp),    dimension(2,2,  ncol,nlay,nflav), intent(out) :: fminor
      !! Interpolation fraction in eta, strat/trop 
    ! -----------------
    ! local
    real(wp), dimension(ncol,nlay) :: ftemp, fpress ! interpolation fraction for temperature, pressure
    real(wp) :: locpress       ! needed to find location in pressure grid
    real(wp) :: ratio_eta_half ! ratio of vmrs of major species that defines eta=0.5
                               ! for given flavor and reference temperature level
    real(wp) :: eta, feta      ! binary_species_parameter, interpolation variable for eta
    real(wp) :: loceta         ! needed to find location in eta grid
    real(wp) :: ftemp_term
    ! -----------------
    ! local indexes
    integer :: icol, ilay, iflav, igases(2), itropo, itemp

    do ilay = 1, nlay
      do icol = 1, ncol
        ! index and factor for temperature interpolation
        jtemp(icol,ilay) = int((tlay(icol,ilay) - (temp_ref_min - temp_ref_delta)) / temp_ref_delta)
        jtemp(icol,ilay) = min(ntemp - 1, max(1, jtemp(icol,ilay))) ! limit the index range
        ftemp(icol,ilay) = (tlay(icol,ilay) - temp_ref(jtemp(icol,ilay))) / temp_ref_delta

        ! index and factor for pressure interpolation
        locpress = 1._wp + (log(play(icol,ilay)) - press_ref_log(1)) / press_ref_log_delta
        jpress(icol,ilay) = min(npres-1, max(1, int(locpress)))
        fpress(icol,ilay) = locpress - float(jpress(icol,ilay))

        ! determine if in lower or upper part of atmosphere
        tropo(icol,ilay) = log(play(icol,ilay)) > press_ref_trop_log
      end do
    end do

    do iflav = 1, nflav
      igases(:) = flavor(:,iflav)
      do ilay = 1, nlay
        do icol = 1, ncol
        ! itropo = 1 lower atmosphere; itropo = 2 upper atmosphere
        itropo = merge(1,2,tropo(icol,ilay))
        ! loop over implemented combinations of major species
          do itemp = 1, 2
            ! compute interpolation fractions needed for lower, then upper reference temperature level
            ! compute binary species parameter (eta) for flavor and temperature and
            !  associated interpolation index and factors
            ratio_eta_half = vmr_ref(itropo,igases(1),(jtemp(icol,ilay)+itemp-1)) / &
                             vmr_ref(itropo,igases(2),(jtemp(icol,ilay)+itemp-1))
            col_mix(itemp,icol,ilay,iflav) = col_gas(icol,ilay,igases(1)) + ratio_eta_half * col_gas(icol,ilay,igases(2))
            ! Keep this commented lines. Fortran does allow for
            ! substantial optimizations and in this merge cases may
            ! happen that all expressions are evaluated and so create
            ! a division by zero. In the if construct this should be
            ! save. Merge is the way to do it in general inside of
            ! loops, but sometimes it may not work.
            !
            ! eta = merge(col_gas(icol,ilay,igases(1)) / col_mix(itemp,icol,ilay,iflav), 0.5_wp, &
            !             col_mix(itemp,icol,ilay,iflav) > 2._wp * tiny(col_mix))
            !
            ! In essence: do not turn it back to merge(...)!
            if (col_mix(itemp,icol,ilay,iflav) > 2._wp * tiny(col_mix)) then
              eta = col_gas(icol,ilay,igases(1)) / col_mix(itemp,icol,ilay,iflav)
            else
              eta = 0.5_wp
            endif
            loceta = eta * float(neta-1)
            jeta(itemp,icol,ilay,iflav) = min(int(loceta)+1, neta-1)
            feta = mod(loceta, 1.0_wp)
            ! compute interpolation fractions needed for minor species
            ! ftemp_term = (1._wp-ftemp(icol,ilay)) for itemp = 1, ftemp(icol,ilay) for itemp=2
            ftemp_term = (real(2-itemp, wp) + real(2*itemp-3, wp) * ftemp(icol,ilay))
            fminor(1,itemp,icol,ilay,iflav) = (1._wp-feta) * ftemp_term
            fminor(2,itemp,icol,ilay,iflav) =        feta  * ftemp_term
            ! compute interpolation fractions needed for major species
            fmajor(1,1,itemp,icol,ilay,iflav) = (1._wp-fpress(icol,ilay)) * fminor(1,itemp,icol,ilay,iflav)
            fmajor(2,1,itemp,icol,ilay,iflav) = (1._wp-fpress(icol,ilay)) * fminor(2,itemp,icol,ilay,iflav)
            fmajor(1,2,itemp,icol,ilay,iflav) =        fpress(icol,ilay)  * fminor(1,itemp,icol,ilay,iflav)
            fmajor(2,2,itemp,icol,ilay,iflav) =        fpress(icol,ilay)  * fminor(2,itemp,icol,ilay,iflav)
          end do ! reference temperatures
        end do ! icol
      end do ! ilay
    end do ! iflav

  end subroutine interpolation
  ! --------------------------------------------------------------------------------------
  !
  !> Compute minor and major species optical depth using pre-computed interpolation coefficients
  !>   (jeta,jtemp,jpress) and weights (fmajor, fminor)
  !
  subroutine compute_tau_absorption(                &
                ncol,nlay,nbnd,ngpt,                &  ! dimensions
                ngas,nflav,neta,npres,ntemp,        &
                nminorlower, nminorklower,          & ! number of minor contributors, total num absorption coeffs
                nminorupper, nminorkupper,          &
                idx_h2o,                            &
                gpoint_flavor,                      &
                band_lims_gpt,                      &
                kmajor,                             &
                kminor_lower,                       &
                kminor_upper,                       &
                minor_limits_gpt_lower,             &
                minor_limits_gpt_upper,             &
                minor_scales_with_density_lower,    &
                minor_scales_with_density_upper,    &
                scale_by_complement_lower,          &
                scale_by_complement_upper,          &
                idx_minor_lower,                    &
                idx_minor_upper,                    &
                idx_minor_scaling_lower,            &
                idx_minor_scaling_upper,            &
                kminor_start_lower,                 &
                kminor_start_upper,                 &
                tropo,                              &
                col_mix,fmajor,fminor,              &
                play,tlay,col_gas,                  &
                jeta,jtemp,jpress,                  &
                tau) bind(C, name="rrtmgp_compute_tau_absorption")
    ! ---------------------
    ! input dimensions
    integer,                                intent(in) :: ncol,nlay,nbnd,ngpt         !! array sizes 
    integer,                                intent(in) :: ngas,nflav,neta,npres,ntemp !! tables sizes 
    integer,                                intent(in) :: nminorlower, nminorklower,nminorupper, nminorkupper
                                                          !! table sizes
    integer,                                intent(in) :: idx_h2o                     !! index of water vapor in col_gas
    ! ---------------------
    ! inputs from object
    integer,     dimension(2,ngpt),                  intent(in) :: gpoint_flavor
      !! major gas flavor (pair) by upper/lower, g-point
    integer,     dimension(2,nbnd),                  intent(in) :: band_lims_gpt
      !! beginning and ending g-point for each band 
    real(wp),    dimension(ntemp,neta,npres+1,ngpt), intent(in) :: kmajor
      !! absorption coefficient table - major gases 
    real(wp),    dimension(ntemp,neta,nminorklower), intent(in) :: kminor_lower
      !! absorption coefficient table - minor gases, lower atmosphere 
    real(wp),    dimension(ntemp,neta,nminorkupper), intent(in) :: kminor_upper
      !! absorption coefficient table - minor gases, upper atmosphere 
    integer,     dimension(2,nminorlower),           intent(in) :: minor_limits_gpt_lower
      !! beginning and ending g-point for each minor gas 
    integer,     dimension(2,nminorupper),           intent(in) :: minor_limits_gpt_upper
    logical(wl), dimension(  nminorlower),           intent(in) :: minor_scales_with_density_lower
      !! generic treatment of minor gases - scales with density (e.g. continuum, collision-induced absorption)?
    logical(wl), dimension(  nminorupper),           intent(in) :: minor_scales_with_density_upper
    logical(wl), dimension(  nminorlower),           intent(in) :: scale_by_complement_lower
      !! generic treatment of minor gases - scale by density (e.g. self-continuum) or complement?  
    logical(wl), dimension(  nminorupper),           intent(in) :: scale_by_complement_upper
    integer,     dimension(  nminorlower),           intent(in) :: idx_minor_lower  
      !! index of each minor gas in col_gas
    integer,     dimension(  nminorupper),           intent(in) :: idx_minor_upper
    integer,     dimension(  nminorlower),           intent(in) :: idx_minor_scaling_lower 
      !! for this minor gas, index of the "scaling gas" in col_gas 
    integer,     dimension(  nminorupper),           intent(in) :: idx_minor_scaling_upper
    integer,     dimension(  nminorlower),           intent(in) :: kminor_start_lower 
      !! starting g-point index in minor gas absorption table
    integer,     dimension(  nminorupper),           intent(in) :: kminor_start_upper
    logical(wl), dimension(ncol,nlay),               intent(in) :: tropo
      !! use upper- or lower-atmospheric tables? 
    ! ---------------------
    ! inputs from profile or parent function
    real(wp), dimension(2,    ncol,nlay,nflav       ), intent(in) :: col_mix
      !! combination of major species's column amounts - computed in interpolation() 
    real(wp), dimension(2,2,2,ncol,nlay,nflav       ), intent(in) :: fmajor
      !! interpolation weights for major gases - computed in interpolation() 
    real(wp), dimension(2,2,  ncol,nlay,nflav       ), intent(in) :: fminor
      !! interpolation weights for minor gases - computed in interpolation() 
    real(wp), dimension(            ncol,nlay       ), intent(in) :: play, tlay 
      !! input temperature and pressure 
    real(wp), dimension(            ncol,nlay,0:ngas), intent(in) :: col_gas
      !! input column gas amount (molecules/cm^2) 
    integer,  dimension(2,    ncol,nlay,nflav       ), intent(in) :: jeta
      !! interpolation indexes in eta - computed in interpolation() 
    integer,  dimension(            ncol,nlay       ), intent(in) :: jtemp
      !! interpolation indexes in temperature - computed in interpolation() 
    integer,  dimension(            ncol,nlay       ), intent(in) :: jpress
      !! interpolation indexes in pressure  - computed in interpolation() 
    ! ---------------------
    ! output - optical depth
    real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau !! aborption optional depth 
    ! ---------------------
    ! Local variables
    !
    logical                    :: top_at_1
    integer, dimension(ncol,2) :: itropo_lower, itropo_upper
    ! ----------------------------------------------------------------

    ! ---------------------
    ! Layer limits of upper, lower atmospheres
    ! ---------------------
    top_at_1 = play(1,1) < play(1, nlay)
    if(top_at_1) then
      itropo_lower(:, 1) = minloc(play, dim=2, mask=tropo)
      itropo_lower(:, 2) = nlay
      itropo_upper(:, 1) = 1
      itropo_upper(:, 2) = maxloc(play, dim=2, mask=(.not. tropo))
    else
      itropo_lower(:, 1) = 1
      itropo_lower(:, 2) = minloc(play, dim=2, mask= tropo)
      itropo_upper(:, 1) = maxloc(play, dim=2, mask=(.not. tropo))
      itropo_upper(:, 2) = nlay
    end if
    ! ---------------------
    ! Major Species
    ! ---------------------
    call gas_optical_depths_major(   &
          ncol,nlay,nbnd,ngpt,       & ! dimensions
          nflav,neta,npres,ntemp,    &
          gpoint_flavor,             &
          band_lims_gpt,             &
          kmajor,                    &
          col_mix,fmajor,            &
          jeta,tropo,jtemp,jpress,   &
          tau)
    ! ---------------------
    ! Minor Species - lower
    ! ---------------------
    call gas_optical_depths_minor(     &
           ncol,nlay,ngpt,             & ! dimensions
           ngas,nflav,ntemp,neta,      &
           nminorlower,nminorklower,   &
           idx_h2o,                    &
           gpoint_flavor(1,:),         &
           kminor_lower,               &
           minor_limits_gpt_lower,     &
           minor_scales_with_density_lower, &
           scale_by_complement_lower,  &
           idx_minor_lower,            &
           idx_minor_scaling_lower,    &
           kminor_start_lower,         &
           play, tlay,                 &
           col_gas,fminor,jeta,        &
           itropo_lower,jtemp,         &
           tau)
    ! ---------------------
    ! Minor Species - upper
    ! ---------------------
    call gas_optical_depths_minor(     &
           ncol,nlay,ngpt,             & ! dimensions
           ngas,nflav,ntemp,neta,      &
           nminorupper,nminorkupper,   &
           idx_h2o,                    &
           gpoint_flavor(2,:),         &
           kminor_upper,               &
           minor_limits_gpt_upper,     &
           minor_scales_with_density_upper, &
           scale_by_complement_upper,  &
           idx_minor_upper,            &
           idx_minor_scaling_upper,    &
           kminor_start_upper,         &
           play, tlay,                 &
           col_gas,fminor,jeta,        &
           itropo_upper,jtemp,         &
           tau)
  end subroutine compute_tau_absorption
  ! --------------------------------------------------------------------------------------

  ! --------------------------------------------------------------------------------------
  !
  ! compute minor species optical depths
  !
  subroutine gas_optical_depths_major(ncol,nlay,nbnd,ngpt,&
                                      nflav,neta,npres,ntemp,      & ! dimensions
                                      gpoint_flavor, band_lims_gpt,   & ! inputs from object
                                      kmajor,                         &
                                      col_mix,fmajor,                 &
                                      jeta,tropo,jtemp,jpress,        & ! local input
                                      tau)
    ! input dimensions
    integer, intent(in) :: ncol, nlay, nbnd, ngpt, nflav,neta,npres,ntemp  ! dimensions

    ! inputs from object
    integer,  dimension(2,ngpt),  intent(in) :: gpoint_flavor
    integer,  dimension(2,nbnd),  intent(in) :: band_lims_gpt ! start and end g-point for each band
    real(wp), dimension(ntemp,neta,npres+1,ngpt), intent(in) :: kmajor

    ! inputs from profile or parent function
    real(wp),    dimension(2,    ncol,nlay,nflav), intent(in) :: col_mix
    real(wp),    dimension(2,2,2,ncol,nlay,nflav), intent(in) :: fmajor
    integer,     dimension(2,    ncol,nlay,nflav), intent(in) :: jeta
    logical(wl), dimension(ncol,nlay), intent(in) :: tropo
    integer,     dimension(ncol,nlay), intent(in) :: jtemp, jpress

    ! outputs
    real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau
    ! -----------------
    ! local variables
    real(wp) :: tau_major(ngpt) ! major species optical depth
    ! local index
    integer :: icol, ilay, iflav, ibnd, itropo
    integer :: gptS, gptE

    ! optical depth calculation for major species
    do ibnd = 1, nbnd
      gptS = band_lims_gpt(1, ibnd)
      gptE = band_lims_gpt(2, ibnd)
      do ilay = 1, nlay
        do icol = 1, ncol
          ! itropo = 1 lower atmosphere; itropo = 2 upper atmosphere
          itropo = merge(1,2,tropo(icol,ilay))
          iflav = gpoint_flavor(itropo, gptS) !eta interpolation depends on band's flavor
          tau_major(gptS:gptE) = &
            ! interpolation in temperature, pressure, and eta
            interpolate3D_byflav(col_mix(:,icol,ilay,iflav),                                     &
                                 fmajor(:,:,:,icol,ilay,iflav), kmajor,                          &
                                 band_lims_gpt(1, ibnd), band_lims_gpt(2, ibnd),                 &
                                 jeta(:,icol,ilay,iflav), jtemp(icol,ilay),jpress(icol,ilay)+itropo)
            tau(icol,ilay,gptS:gptE) = tau(icol,ilay,gptS:gptE) + tau_major(gptS:gptE)
        end do
      end do
    end do

  end subroutine gas_optical_depths_major

  ! ----------------------------------------------------------
  !
  ! compute minor species optical depths
  !
  subroutine gas_optical_depths_minor(ncol,nlay,ngpt,        &
                                      ngas,nflav,ntemp,neta, &
                                      nminor,nminork,        &
                                      idx_h2o,               &
                                      gpt_flv,               &
                                      kminor,                &
                                      minor_limits_gpt,      &
                                      minor_scales_with_density,    &
                                      scale_by_complement,   &
                                      idx_minor, idx_minor_scaling, &
                                      kminor_start,          &
                                      play, tlay,            &
                                      col_gas,fminor,jeta,   &
                                      layer_limits,jtemp,    &
                                      tau)
    integer,                                     intent(in   ) :: ncol,nlay,ngpt
    integer,                                     intent(in   ) :: ngas,nflav
    integer,                                     intent(in   ) :: ntemp,neta,nminor,nminork
    integer,                                     intent(in   ) :: idx_h2o
    integer,     dimension(ngpt),                intent(in   ) :: gpt_flv
    real(wp),    dimension(ntemp,neta,nminork),  intent(in   ) :: kminor
    integer,     dimension(2,nminor),            intent(in   ) :: minor_limits_gpt
    logical(wl), dimension(  nminor),            intent(in   ) :: minor_scales_with_density
    logical(wl), dimension(  nminor),            intent(in   ) :: scale_by_complement
    integer,     dimension(  nminor),            intent(in   ) :: kminor_start
    integer,     dimension(  nminor),            intent(in   ) :: idx_minor, idx_minor_scaling
    real(wp),    dimension(ncol,nlay),           intent(in   ) :: play, tlay
    real(wp),    dimension(ncol,nlay,0:ngas),    intent(in   ) :: col_gas
    real(wp),    dimension(2,2,ncol,nlay,nflav), intent(in   ) :: fminor
    integer,     dimension(2,  ncol,nlay,nflav), intent(in   ) :: jeta
    integer,     dimension(ncol, 2),             intent(in   ) :: layer_limits
    integer,     dimension(ncol,nlay),           intent(in   ) :: jtemp
    real(wp),    dimension(ncol,nlay,ngpt),      intent(inout) :: tau
    ! -----------------
    ! local variables
    real(wp), parameter :: PaTohPa = 0.01_wp
    real(wp) :: vmr_fact, dry_fact             ! conversion from column abundance to dry vol. mixing ratio;
    real(wp) :: scaling                        ! optical depth
    integer  :: icol, ilay, iflav, imnr
    integer  :: gptS, gptE
    real(wp), dimension(ngpt) :: tau_minor
    ! -----------------
    !
    ! Guard against layer limits being 0 -- that means don't do anything i.e. there are no
    !   layers with pressures in the upper or lower atmosphere respectively
    ! First check skips the routine entirely if all columns are out of bounds...
    !

    if(any(layer_limits(:,1) > 0)) then
      do imnr = 1, size(scale_by_complement,dim=1) ! loop over minor absorbers in each band
        do icol = 1, ncol
          !
          ! This check skips individual columns with no pressures in range
          !
          if(layer_limits(icol,1) > 0) then
            do ilay = layer_limits(icol,1), layer_limits(icol,2)
              !
              ! Scaling of minor gas absortion coefficient begins with column amount of minor gas
              !
              scaling = col_gas(icol,ilay,idx_minor(imnr))
              !
              ! Density scaling (e.g. for h2o continuum, collision-induced absorption)
              !
              if (minor_scales_with_density(imnr)) then
                !
                ! NOTE: P needed in hPa to properly handle density scaling.
                !
                scaling = scaling * (PaTohPa*play(icol,ilay)/tlay(icol,ilay))
                if(idx_minor_scaling(imnr) > 0) then  ! there is a second gas that affects this gas's absorption
                  vmr_fact = 1._wp / col_gas(icol,ilay,0)
                  dry_fact = 1._wp / (1._wp + col_gas(icol,ilay,idx_h2o) * vmr_fact)
                  ! scale by density of special gas
                  if (scale_by_complement(imnr)) then ! scale by densities of all gases but the special one
                    scaling = scaling * (1._wp - col_gas(icol,ilay,idx_minor_scaling(imnr)) * vmr_fact * dry_fact)
                  else
                    scaling = scaling *         (col_gas(icol,ilay,idx_minor_scaling(imnr)) * vmr_fact * dry_fact)
                  endif
                endif
              endif
              !
              ! Interpolation of absorption coefficient and calculation of optical depth
              !
              ! Which gpoint range does this minor gas affect?
              gptS = minor_limits_gpt(1,imnr)
              gptE = minor_limits_gpt(2,imnr)
              iflav = gpt_flv(gptS)
              tau_minor(gptS:gptE) = scaling *                   &
                                      interpolate2D_byflav(fminor(:,:,icol,ilay,iflav), &
                                                           kminor, &
                                                           kminor_start(imnr), kminor_start(imnr)+(gptE-gptS), &
                                                           jeta(:,icol,ilay,iflav), jtemp(icol,ilay))
              tau(icol,ilay,gptS:gptE) = tau(icol,ilay,gptS:gptE) + tau_minor(gptS:gptE)
            enddo
          end if
        enddo
      enddo
    end if

  end subroutine gas_optical_depths_minor
  ! ----------------------------------------------------------
  !
  ! compute Rayleigh scattering optical depths
  !
  subroutine compute_tau_rayleigh(ncol,nlay,nbnd,ngpt,         &
                                  ngas,nflav,neta,npres,ntemp, &
                                  gpoint_flavor,band_lims_gpt, &
                                  krayl,                       &
                                  idx_h2o, col_dry,col_gas,    &
                                  fminor,jeta,tropo,jtemp,     &
                                  tau_rayleigh) bind(C, name="rrtmgp_compute_tau_rayleigh")
    integer,                                     intent(in ) :: ncol,nlay,nbnd,ngpt
      !! input dimensions 
    integer,                                     intent(in ) :: ngas,nflav,neta,npres,ntemp
      !! table dimensions 
    integer,     dimension(2,ngpt),              intent(in ) :: gpoint_flavor 
      !! major gas flavor (pair) by upper/lower, g-point
    integer,     dimension(2,nbnd),              intent(in ) :: band_lims_gpt
      !! start and end g-point for each band
    real(wp),    dimension(ntemp,neta,ngpt,2),   intent(in ) :: krayl
      !! Rayleigh scattering coefficients 
    integer,                                     intent(in ) :: idx_h2o
      !! index of water vapor in col_gas
    real(wp),    dimension(ncol,nlay),           intent(in ) :: col_dry
      !! column amount of dry air 
    real(wp),    dimension(ncol,nlay,0:ngas),    intent(in ) :: col_gas
      !! input column gas amount  (molecules/cm^2)
    real(wp),    dimension(2,2,ncol,nlay,nflav), intent(in ) :: fminor
      !! interpolation weights for major gases - computed in interpolation() 
    integer,     dimension(2,  ncol,nlay,nflav), intent(in ) :: jeta
      !! interpolation indexes in eta - computed in interpolation() 
    logical(wl), dimension(ncol,nlay),           intent(in ) :: tropo
      !! use upper- or lower-atmospheric tables? 
    integer,     dimension(ncol,nlay),           intent(in ) :: jtemp
      !! interpolation indexes in temperature - computed in interpolation() 
    ! outputs
    real(wp),    dimension(ncol,nlay,ngpt),      intent(out) :: tau_rayleigh
      !! Rayleigh optical depth 
    ! -----------------
    ! local variables
    real(wp) :: k(ngpt) ! rayleigh scattering coefficient
    integer  :: icol, ilay, iflav, ibnd, gptS, gptE
    integer  :: itropo
    ! -----------------

    do ibnd = 1, nbnd
      gptS = band_lims_gpt(1, ibnd)
      gptE = band_lims_gpt(2, ibnd)
      do ilay = 1, nlay
        do icol = 1, ncol
          itropo = merge(1,2,tropo(icol,ilay)) ! itropo = 1 lower atmosphere;itropo = 2 upper atmosphere
          iflav = gpoint_flavor(itropo, gptS) !eta interpolation depends on band's flavor
          k(gptS:gptE) = interpolate2D_byflav(fminor(:,:,icol,ilay,iflav), &
                                              krayl(:,:,:,itropo),      &
                                              gptS, gptE, jeta(:,icol,ilay,iflav), jtemp(icol,ilay))
          tau_rayleigh(icol,ilay,gptS:gptE) = k(gptS:gptE) * &
                                              (col_gas(icol,ilay,idx_h2o)+col_dry(icol,ilay))
        end do
      end do
    end do

  end subroutine compute_tau_rayleigh

  ! ----------------------------------------------------------
  subroutine compute_Planck_source(                        &
                    ncol, nlay, nbnd, ngpt,                &
                    nflav, neta, npres, ntemp, nPlanckTemp,&
                    tlay, tlev, tsfc, sfc_lay,             &
                    fmajor, jeta, tropo, jtemp, jpress,    &
                    gpoint_bands, band_lims_gpt,           &
                    pfracin, temp_ref_min, totplnk_delta, totplnk, gpoint_flavor, &
                    sfc_src, lay_src, lev_src_inc, lev_src_dec, sfc_source_Jac) bind(C, name="rrtmgp_compute_Planck_source")
    integer,                                    intent(in) :: ncol, nlay, nbnd, ngpt
      !! input dimensions 
    integer,                                    intent(in) :: nflav, neta, npres, ntemp, nPlanckTemp
      !! table dimensions 
    real(wp),    dimension(ncol,nlay  ),        intent(in) :: tlay !! temperature at layer centers (K)
    real(wp),    dimension(ncol,nlay+1),        intent(in) :: tlev !! temperature at interfaces (K)
    real(wp),    dimension(ncol       ),        intent(in) :: tsfc !! surface temperture 
    integer,                                    intent(in) :: sfc_lay !! index into surface layer 
    ! Interpolation variables
    real(wp),    dimension(2,2,2,ncol,nlay,nflav), intent(in) :: fmajor
      !! interpolation weights for major gases - computed in interpolation() 
    integer,     dimension(2,    ncol,nlay,nflav), intent(in) :: jeta
      !! interpolation indexes in eta - computed in interpolation() 
    logical(wl), dimension(            ncol,nlay), intent(in) :: tropo
      !! use upper- or lower-atmospheric tables? 
    integer,     dimension(            ncol,nlay), intent(in) :: jtemp, jpress
      !! interpolation indexes in temperature and pressure - computed in interpolation() 
    ! Table-specific
    integer, dimension(ngpt),                     intent(in) :: gpoint_bands  !! band to which each g-point belongs
    integer, dimension(2, nbnd),                  intent(in) :: band_lims_gpt !! start and end g-point for each band
    real(wp),                                     intent(in) :: temp_ref_min, totplnk_delta !! interpolation constants
    real(wp), dimension(ntemp,neta,npres+1,ngpt), intent(in) :: pfracin       !! Fraction of the Planck function in each g-point
    real(wp), dimension(nPlanckTemp,nbnd),        intent(in) :: totplnk       !! Total Planck function by band at each temperature 
    integer,  dimension(2,ngpt),                  intent(in) :: gpoint_flavor !! major gas flavor (pair) by upper/lower, g-point

    real(wp), dimension(ncol,     ngpt), intent(out) :: sfc_src  !! Planck emssion from the surface 
    real(wp), dimension(ncol,nlay,ngpt), intent(out) :: lay_src  !! Planck emssion from layer centers
    real(wp), dimension(ncol,nlay,ngpt), intent(out) :: lev_src_inc, lev_src_dec
      !! Planck emission at layer boundaries, using spectral mapping in the direction of propagation 
    real(wp), dimension(ncol,     ngpt), intent(out) :: sfc_source_Jac 
      !! Jacobian (derivative) of the surface Planck source with respect to surface temperature 
    ! -----------------
    ! local
    real(wp), parameter                             :: delta_Tsurf = 1.0_wp

    integer  :: ilay, icol, igpt, ibnd, itropo, iflav
    integer  :: gptS, gptE
    real(wp), dimension(2), parameter :: one = [1._wp, 1._wp]
    real(wp) :: pfrac          (ncol,nlay  ,ngpt)
    real(wp) :: planck_function(ncol,nlay+1,nbnd)
    ! -----------------

    ! Calculation of fraction of band's Planck irradiance associated with each g-point
    do ibnd = 1, nbnd
      gptS = band_lims_gpt(1, ibnd)
      gptE = band_lims_gpt(2, ibnd)
      do ilay = 1, nlay
        do icol = 1, ncol
          ! itropo = 1 lower atmosphere; itropo = 2 upper atmosphere
          itropo = merge(1,2,tropo(icol,ilay))
          iflav = gpoint_flavor(itropo, gptS) !eta interpolation depends on band's flavor
          pfrac(icol,ilay,gptS:gptE) = &
            ! interpolation in temperature, pressure, and eta
            interpolate3D_byflav(one, fmajor(:,:,:,icol,ilay,iflav), pfracin, &
                          band_lims_gpt(1, ibnd), band_lims_gpt(2, ibnd),                 &
                          jeta(:,icol,ilay,iflav), jtemp(icol,ilay),jpress(icol,ilay)+itropo)
        end do ! column
      end do   ! layer
    end do     ! band

    !
    ! Planck function by band for the surface
    ! Compute surface source irradiance for g-point, equals band irradiance x fraction for g-point
    !
    do icol = 1, ncol
      planck_function(icol,1,1:nbnd) = interpolate1D(tsfc(icol), temp_ref_min, totplnk_delta, totplnk)
      planck_function(icol,2,1:nbnd) = interpolate1D(tsfc(icol) + delta_Tsurf, temp_ref_min, totplnk_delta, totplnk)
      !
      ! Map to g-points
      !
      do ibnd = 1, nbnd
        gptS = band_lims_gpt(1, ibnd)
        gptE = band_lims_gpt(2, ibnd)
        do igpt = gptS, gptE
            sfc_src(icol,igpt) = pfrac(icol,sfc_lay,igpt) * planck_function(icol,1,ibnd)
            sfc_source_Jac(icol, igpt) = pfrac(icol,sfc_lay,igpt) * &
                                (planck_function(icol, 2, ibnd) - planck_function(icol,1,ibnd))
        end do
      end do
    end do !icol

    do ilay = 1, nlay
      do icol = 1, ncol
        ! Compute layer source irradiance for g-point, equals band irradiance x fraction for g-point
        planck_function(icol,ilay,1:nbnd) = interpolate1D(tlay(icol,ilay), temp_ref_min, totplnk_delta, totplnk)
      end do
    end do

    !
    ! Map to g-points
    !
    do ibnd = 1, nbnd
      gptS = band_lims_gpt(1, ibnd)
      gptE = band_lims_gpt(2, ibnd)
      do igpt = gptS, gptE
        do ilay = 1, nlay
          do icol = 1, ncol
            lay_src(icol,ilay,igpt) = pfrac(icol,ilay,igpt) * planck_function(icol,ilay,ibnd)
          end do
        end do
      end do
    end do

    ! compute level source irradiances for each g-point, one each for upward and downward paths
    do ilay = 1, nlay
      do icol = 1, ncol
      planck_function(icol,     1,1:nbnd) = interpolate1D(tlev(icol,     1),temp_ref_min, totplnk_delta, totplnk)
      planck_function(icol,ilay+1,1:nbnd) = interpolate1D(tlev(icol,ilay+1),temp_ref_min, totplnk_delta, totplnk)
      end do
    end do

    !
    ! Map to g-points
    !
    do ibnd = 1, nbnd
      gptS = band_lims_gpt(1, ibnd)
      gptE = band_lims_gpt(2, ibnd)
      do igpt = gptS, gptE
        do ilay = 1, nlay
          do icol = 1, ncol
            lev_src_inc(icol,ilay,igpt) = pfrac(icol,ilay,igpt) *planck_function(icol,ilay+1,ibnd)
            lev_src_dec(icol,ilay,igpt) = pfrac(icol,ilay,igpt) *planck_function(icol,ilay  ,ibnd)
          end do
        end do
      end do
    end do

  end subroutine compute_Planck_source
  ! ----------------------------------------------------------
  !
  ! One dimensional interpolation -- return all values along second table dimension
  !
  pure function interpolate1D(val, offset, delta, table) result(res)
    ! input
    real(wp), intent(in) :: val,    & ! axis value at which to evaluate table
                            offset, & ! minimum of table axis
                            delta     ! step size of table axis
    real(wp), dimension(:,:), &
              intent(in) :: table ! dimensions (axis, values)
    ! output
    real(wp), dimension(size(table,dim=2)) :: res

    ! local
    real(wp) :: val0 ! fraction index adjusted by offset and delta
    integer :: index ! index term
    real(wp) :: frac ! fractional term
    ! -------------------------------------
    val0 = (val - offset) / delta
    frac = val0 - int(val0) ! get fractional part
    index = min(size(table,dim=1)-1, max(1, int(val0)+1)) ! limit the index range
    res(:) = table(index,:) + frac * (table(index+1,:) - table(index,:))
  end function interpolate1D
  ! ----------------------------------------------------------
  !   This function returns a range of values from a subset (in gpoint) of the k table
  !
  pure function interpolate2D_byflav(fminor, k, gptS, gptE, jeta, jtemp) result(res)
    real(wp), dimension(2,2), intent(in) :: fminor ! interpolation fractions for minor species
                                       ! index(1) : reference eta level (temperature dependent)
                                       ! index(2) : reference temperature level
    real(wp), dimension(:,:,:), intent(in) :: k ! (g-point, eta, temp)
    integer,                    intent(in) :: gptS, gptE, jtemp ! interpolation index for temperature
    integer, dimension(2),      intent(in) :: jeta ! interpolation index for binary species parameter (eta)
    real(wp), dimension(gptE-gptS+1)       :: res ! the result

    ! Local variable
    integer :: igpt
    ! each code block is for a different reference temperature

    do igpt = 1, gptE-gptS+1
      res(igpt) = fminor(1,1) * k(jtemp  , jeta(1)  , gptS+igpt-1) + &
                  fminor(2,1) * k(jtemp  , jeta(1)+1, gptS+igpt-1) + &
                  fminor(1,2) * k(jtemp+1, jeta(2)  , gptS+igpt-1) + &
                  fminor(2,2) * k(jtemp+1, jeta(2)+1, gptS+igpt-1)
    end do

  end function interpolate2D_byflav
  ! ----------------------------------------------------------
  pure function interpolate3D_byflav(scaling, fmajor, k, gptS, gptE, jeta, jtemp, jpress) result(res)
    real(wp), dimension(2),     intent(in) :: scaling
    real(wp), dimension(2,2,2), intent(in) :: fmajor ! interpolation fractions for major species
                                                     ! index(1) : reference eta level (temperature dependent)
                                                     ! index(2) : reference pressure level
                                                     ! index(3) : reference temperature level
    real(wp), dimension(:,:,:,:),intent(in) :: k ! (temp,eta,press,gpt)
    integer,                     intent(in) :: gptS, gptE
    integer, dimension(2),       intent(in) :: jeta ! interpolation index for binary species parameter (eta)
    integer,                     intent(in) :: jtemp ! interpolation index for temperature
    integer,                     intent(in) :: jpress ! interpolation index for pressure
    real(wp), dimension(gptS:gptE)          :: res ! the result

    ! Local variable
    integer :: igpt
    ! each code block is for a different reference temperature
    do igpt = gptS, gptE
      res(igpt) =  &
        scaling(1) * &
        ( fmajor(1,1,1) * k(jtemp, jeta(1)  , jpress-1, igpt) + &
          fmajor(2,1,1) * k(jtemp, jeta(1)+1, jpress-1, igpt) + &
          fmajor(1,2,1) * k(jtemp, jeta(1)  , jpress  , igpt) + &
          fmajor(2,2,1) * k(jtemp, jeta(1)+1, jpress  , igpt) ) + &
        scaling(2) * &
        ( fmajor(1,1,2) * k(jtemp+1, jeta(2)  , jpress-1, igpt) + &
          fmajor(2,1,2) * k(jtemp+1, jeta(2)+1, jpress-1, igpt) + &
          fmajor(1,2,2) * k(jtemp+1, jeta(2)  , jpress  , igpt) + &
          fmajor(2,2,2) * k(jtemp+1, jeta(2)+1, jpress  , igpt) )
    end do
  end function interpolate3D_byflav

end module mo_gas_optics_rrtmgp_kernels
